
function [Q,T]=Schurl(W);

%     This procedure produces the 
%     Schur decomposition of the matrix W
%     That is  Q'*W*Q=T
%
%     Where T is LOWER triangular  (Note lower and not upper)
%     and the diagonal entries of T are the ordered
%     eigenvalues of  W with T[i,i] >= T[j,j], for i < j.



if (rows(W) == 1); 
 T=W; Q=1;
else

%  Step 1:  Find Real Schur Decomposition  %
% description of MATLAB's schur function, which
%       does not find the desired schur decomposition
%
% [U,T] = SCHUR(X) produces a Schur matrix T and a unitary
%       matrix U so that X = U*T*U' and U'*U = EYE(U).
%       By itself, SCHUR(X) returns T.
%       If X is complex, the Complex Schur Form is returned in matrix T.
%       The Complex Schur Form is upper triangular with the eigenvalues
%       of X on the diagonal.
%       If X is real, the Real Schur Form is returned.  The Real Schur
%       Form has the real eigenvalues on the diagonal and the complex
%       eigenvalues in 2-by-2 blocks on the diagonal.
%       See RSF2CSF to convert from Real to Complex Schur form.

 
[Q,T]=schur(W);

%  Step 2:  Convert to Strictly Upper Triangular, possibly Complex Form %
%
%  description of MATLAB's real complex schur form conversion
%
%  RSF2CSF Real block diagonal form to complex diagonal form.
%       [U,T] = RSF2CSF(U,T) converts a real, upper quasi-triangular
%       Schur form to a complex, upper triangular Schur form.

[Q,T]=rsf2csf(Q,T);

% Step 3:  Order Diagonal Elements of T
% from Smallest to Largest 

i=2;

while i <= rows(T);
  if abs(T(i-1,i-1)) > abs(T(i,i));
     [T,Q]=tqswitch(T,Q,i-1);
   if i > 2; i=i-2; end;
  end;
i=i+1;
end;

% Step 4:  Reorder Elements So that T is Upper Triangular%
Q1=Q;
i=1;

while i <= cols(Q);
 Q(:,i)=Q1(:,cols(Q1)+1-i);
 i=i+1; 
end;

T=rev(rev(T)')';

end


